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Abstract 

The aim of this paper is to provide a fast and efficient procedure for (real-time) 
target identification in imaging based on matching on a dictionary of prccomputcd 
generalized polarization tensors (GPTs). The approach is based on some important 
properties of the GPTs and new invariants. A new shape representation is given and 
numerically tested in the presence of measurement noise. The stability and resolution 
of the proposed identification algorithm is numerically quantified. 
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1 Introduction 

With each domain and material parameter, an infinite number of tensors, called the Gener- 
alized Polarization Tensors (GPTs), is associated. The concept of GPTs was introduced in 
[8, 6]. The GPTs contain significant information on the shape of the domain [9]. It occurs 
in several interesting contexts, in particular, in low- frequency scattering [17, 6], asymptotic 
models of dilute composites (see [21] and [13]), in invisibility cloaking in the quasi-static 
regime [10] and in potential theory related to certain questions arising in hydrodynamics 
[22]. 

Another important use of this concept is for imaging diametrically small inclusions from 
boundary measurements. In fact, the GPTs are the basic building blocks for the asymptotic 
expansions of the boundary voltage perturbations due to the presence of small conductivity 
inclusions inside a conductor [18, 16, 8]. Based on this expansion, efficient algorithms to 
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determine the location and some geometric features of the inclusions were proposed. We 
refer to [6, 7] and the references therein for recent developments of this theory. 

In [11], a recursive optimal control scheme to recover fine shape details of a given domain 
using GPTs is proposed. In [4], it is shown that high-frequency oscillations of the boundary 
of a domain are only contained in its high-order GPTs. Moreover, by developing a level set 
version of the recursive optimization scheme, it is also shown that the GPTs can capture the 
topology of the domain. An efficient algorithm for computing the GPTs has been presented 
in [15]. 

The aim of this paper is to show that the GPTs can be used for target identification from 
imaging data. In fact, the GPTs can be accurately obtained from multistatic measurements 
by solving a linear system. Based on this, we design a fast algorithm which identifies a 
target using a dictionary of precomputed GPTs data. We first provide a stability analysis 
for the reconstruction of the GPTs in the presence of measurement noise which quantifies 
the ill-posedness of the imaging problem. Then, suppose that we have a dictionary which 
is a collection of standard shapes (for example alphabetic letters or flowers). Our aim is to 
identify from imaging data a shape which is obtained from one element of the dictionary after 
some rotation, scaling and translation. We design a dictionary matching procedure which 
operates directly in the GPTs data. Our procedure is based on some important properties 
of the GPTs and new invariants. We test the robustness of our procedure with respect 
to a measurement noise in the imaging data. Our approach is quite natural since it uses 
geometric quantities obtained from the imaging data by simply inverting a linear system. 
Moreover, there is an infinite number of invariants associated with the GPTs. Furthermore, 
for a given dictionary, the GPT-based representation may lead to better distinguishibility 
between the dictionary elements. 

Over the last decades, a considerable amount of work has been devoted to nonlinear 
optimization techniques for solving the imaging problem; see, for instance, [19, 23, 25] and 
the references therein. More recently, new regularized optimal control formulations for 
target imaging have been proposed in [1, 3]. As far as we know, our approach in this paper 
provides for the first time an alternative approach to solving the full inverse problem for 
target identification and characterization. It opens a way for real-time target identification 
and tracking algorithms in wave imaging. 

The paper is organized as follows. In section 2, we introduce a particular linear combi- 
nation of the GPTs to obtain what we call the contracted GPTs (CGPTs) [10]. In Section 
3, we investigate the reconstruction of contracted GPTs, defined in (2.14)-(2.17) below, 
from the multistatic response matrix of a conductivity problem. We also consider the effect 
of the presence of measurement noise in the MSR on the reconstruction of the CGPTs. 
Given a signal-to-noise ratio, we determine the statistical stability in the reconstruction 
of the CGPTs, and show that such inverse problem is exponentially unstable. This is the 
well-known ill-posedness of the inverse conductivity problem. In section 4 it is shown that 
the CGPTs have some nice properties, such as simple rotation and translation formulas, 
simple relation with shape symmetry, etc. More importantly, we derive new invariants for 
the CGPTs. One of the matching algorithms presented in section 5 is based on those invari- 
ants. Section 6 presents a variety of numerical results for the target identification problem 
and shows the viability of the proposed procedure. 
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2 Structure of the Multistatic Response Matrix 



The first part of this paper is to reconstruct CGPTs from the multistatic response (MSR) 
matrix, which measures the change in potential field due to a conductivity inclusion. In this 
section, we present the mathematical model for MSR and write it in terms of the CGPTs 
associated to the conductivity inclusion. 

We consider a two dimensional conductivity medium with uniform conductivity equal to 
one, except in an inclusion where the conductivity is k > 1; we denote by A the contrast of 
this inclusion, that is, A = (k + 1)/(2(k — 1)). Let D = z + 5B = {x = z + 5y\y£ B} model 
the conductivity inclusion. Here, B is some C 2 and bounded domain in ]R 2 whose typical 
length scale is of order one; z is a point in M 2 and is taken here to be an estimation of the 
location of the inclusion; 5 is the typical length scale of the inclusion. We refer to [14, 6] 
for efficient location search algorithms and to [2] for correcting the effect of measurement 
noise on the localization procedure. 

The MSR matrix is constructed as follows. Let {x r }^ 1 and {a^}^ model a set of 
electric potential point detectors and electric point sources. We assume in this paper that 
the two sets of locations coincide and N r = N s = N. The MSR matrix V is an iV-by- 
N matrix whose rs-element is the difference of electric potentials with and without the 
conductivity inclusions: 

V rs = u s {x r ) - T s (x r ), r,s = l,...,N. (2.1) 

Here, T s (x) = T(x — x s ) and T(x) = ^ log \x\ is the fundamental solution of the Laplace 
equation in K 2 , and u s (x) is the solution to the transmission problem 

f V • (1 + (k - 1) X d)Vu s (x) = 6 Xs (x), x e R 2 \dD, 

u s( x )\ + = u s(x)\_, x G dD, 

v x ■ (Vii s )| + = nv x ■ (V« s )|_, x G dD, 

u s (x) — T s (x) = \x — x s | — y oo. 



(2.2) 



In the second and third equations above, the notation </>| ± (x) denotes the limit lim t |o (/>(x± 
tu x ), where x G dD and v x is the outward unit normal of dD at x. 

2.1 The asymptotic expansion of the perturbed potential field 

As modeled above, the MSR matrix characterizes the perturbed potential field u s (x r ) — 
T s (x r ). In this section we recall, from [6], the asymptotic expansion of this perturbation 
and some key notions along the way. 

Let So be the single layer potential associated with D, that is, 

S D [4>]{x):= f F(x-y)<j>(y)ds(y), x G M 2 , (2.3) 

J 3D 

and let K-d '■ L 2 (dD) — > L 2 (dD) denote the Poincare- Neumann operator 

K D [<l>]{x) := -!- / {y ~ X '^ ) ( f>(y)ds(y), x G dD. (2.4) 
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Here, (, } denotes the scalar product in M? and u y is the unit normal vector along the 
boundary at y. It is well known that the single layer potential 5d[0] is a harmonic function 
satisfying 5d[^]|_ = <5dM| + and the jump condition 



d_ 

du 



S D [ 



(2.5) 



where KL* D is the adjoint operator of Kd and it has a similar expression as (2.4) with the 
numerator of the integrand replaced by (x — y, u x ). Using (2.5), we verify that T s (x)+Sd[4's] 
with (j) s G L 2 (dD) solving 

8T 

(\I-)C* D )[tj, s } = ?-± 



du 



3D 



(2.6) 



is a solution to the transmission problem (2.2). In fact, this solution is unique and we 
conclude that 



u s (x)-r s ( x ) = s D [<p s }= [ v{x - y )(\i - ie D y l 

J 3D 



1 


'dT s 






du 


dD 



(y)ds(y). 



(2.7) 



To verify the formal derivation above, we refer the reader to Section 2.4 of [6]. 

We assume that the inclusion D and the point z is away from the sources. As a result, 
the functions T(x r — y) and T s (y) are smooth for y & D, and the perturbed field (2.7) is 
well defined. For y € dD and z away from x, the K-th order Taylor expansion formula with 
remainder ex states 

r(x -y)=T(x-z-{y-z))=Y j tlt[ d « T{x - z ){y- z) a + e K . (2.8) 



Throughout this section, we use Greek letters to denote double indices: a = (01,02) €= N 2 , 
a\ = ai!«2! and \a\ = a\ + 02- Substitution of this expansion into (2.7) yields the following 
expansion of V rs plus an error term denoted by E rs : 



A 



Vr, 



2 ^r^ 1 ^ " z)Q aP {z)dPY{z - x s ) + E rs , 



l«l,l/9|=l 



with 



3D 



(y - z) a (XI - K* D y l 



d_ 

du 



{■-4 



(y)ds(y). 



The zeroth order term with f3 = vanishes because the differentiation d/du; the zeroth 
order term corresponding to a = vanishes because (XI — K,* D )~ l maps a zero mean value 
function on dD to another zero mean value function. 

For a generic conductivity inclusion D with the contrast factor A, the GPT of order a/3 
associated with the inclusion is defined by 



M a p{\D):= [ y P(\I-1C* D r l [^-y a }ds{y). 

JdD 



(2.9) 
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Using the change of variable y — z i— > y, the integral term Q a p{z) inside the expansion 
of V rs above can be written as 

Q a p(z)= [ r^i-JChrH^-y^dsiv), (2-io) 

Jd(5B) OU 

which is independant of z. Moreover, by the definition of GPT, this term is Mp a {\,5B). 
As a result, we have 

K 1 

V rs = Y -^d a T(z-x s )M a ^X,5B)d' 3 T(z-x r ) + E rs , (2.11) 

where E rs is the truncation error resulted from the finite expansion. Note also that we have 
switched the indices a and j3. 

The MSR matrix V consisting of u s (x r ) — T s (x r ) depends only on the inclusion (A, D). 
However, the GPTs involved in the representation (2.11) depend on the (non-unique) char- 
acterization (z, 5B) of D. We note that the remainder ex and the truncation error E rs can 
be evaluated; see Appendix A.l. Moreover, since the sensors and the receivers coincide, the 
MSR matrix is symmetric; see (A. 2). 

2.2 Expansion for MSR using contracted GPT 

In this section, we further simplify the expression of MSR using the notion of contracted 
GPT (CGPT), which has been introduced in [10]. Using CGPT, we can write the MSR ma- 
trix V as a product of a CGPT matrix with coefficient matrices, which is a very convenient 
form for inversion. 

Let P m (x) be the complex valued polynomial 

P m ( x ) = (a* + ix 2 ) m := E < xC * + * E b f X ^ ( 2 - 12 ) 

|a|=m |/3|=m 

Using polar coordinate x = re t6 , the above coefficients a™ and bT 1 can also be characterized 
by 

E a™x a =r m cosm6, and Y b^x 13 = r m smm6. (2.13) 

|a|=m |o|=m 

For a generic conductivity inclusion D with contrast A, the associated GPT M a p(\,D) is 
defined as in (2.9). The associated CGPT is the following combination of GPTs using the 
coefficients in (2.12): 

M -n =EE <^M^, (2.14) 

|«|=m |/3|=n 

M ^n = E E ( 2 - 15 ) 

o|=m |/3|=n 

M -n = E E b >P M *f>> ( 2 - 16 ) 
a|=m |/3|=n 

M ^n = E E KW«P- ( 2 - 17 ) 
|a|=m |/3|=n 
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Using the complex coordinate x = r x e x , we have (see Appendix A. 2) that 



^ {—d Q T(x) 



a l 



2ir a 



cos a 



\a\ 

r x 



sin \a\t 



(2.18) 



Recall that {x r }^ =l and {x< j }^ 1 denote the locations of the receivers and electric sources. 
Define R r and 9 r so that the complex representation of x r — z is R T e ldr with z being the 
location of the target. Similarly define R s and 6 S . Substituting formula (2.18) into the 
expression (2.11) of the MSR, we get 



K 

v rs = y, 

|a|=l,|/3|: 



M i i 
a a cos a 



6^' sin \a\9 s 



27r|a|i?i a| 



a ] S 1 cos \0\0r + sin |/3|0 r 



K 



1~L,n = l y> 



/ cos mOe sinm^. 



2irmR" 1 2nmR r : 



M cc 

mn 
M sc 



M c 

n 

M s 



cos np, 
sin n# T 



27r|^|i? ; 
1 



r 



+ E r 



2-khRJ} 



+E r 



Here, the short-hand notations M mn , and A sm 



A sm M mn (Arn)' 

(2.19) 

represent the two-by-two and one-by-two 
matrices respectively, and (A rn )* is the transpose. As m, n run from one to K, which is the 
truncation order of CGPT, and r, s run from one to N, which is the number of receivers 
(sources), these matrices build up the 2K x 2K CGPT block matrix M and the N x 2K 
coefficient matrix A as follows: 



M 



M 21 



M12 
M 22 



\M K1 M K2 



M 1K \ 
M 2K 

M KK J 



A 2 i 



A12 

A 22 



\Ajvi A 



N2 



AnkJ 



(2.20) 



Using these notations, the MSR matrix V can be written as 



V = AMA* + E, 



(2.21) 



where A* denotes the transpose of A and the matrix E = (E rs ) represents the truncation 
error. We precise again that the CGPT above is for the "shifted" inclusion 5B. We note 
also that the dimension of V depends on the number of sources/receivers but does not 
depend on the expansion order K in (2.11). 

Due to the symmetry of harmonic combination of GPTs [7], the matrix M is symmetric. 
Since V is symmetric as shown in (A. 2), the truncation error E is also symmetric. 



3 Reconstruction of CGPTs and Stability Analysis 

The first step in the target identification procedure is to reconstruct CGPTs from the MSR 
matrix V, which has expression (2.21). Define the linear operator L : ]R2^x2^ _^ R NxN by 

L(M) := AMA'. (3.1) 
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We reconstruct CGPTs as the least squares solution of the above linear system, i.e., 

M est = min ||V-L(M tcst )|| F , (3.2) 



mm 

M tcst _Lker(L) 



where ker(L) denotes the kernel of L and || • \\p denotes the Frobenius norm of matrices 
[20]. In general we take N large enough so that 2K < N. When A has full rank 2K, L is 
rank preserving and ker(L) is trivial; in that case, the admissible set above can be replaced 

by R^x2K and 



M 



(A A) A VA(A A) 



From the structure of the matrix A in (2.20) and the expression of the MSR matrix, 
we observe that the contribution of a CGPT decays as its order grows. Consequently, one 
does not expect the inverse procedure to be stable for higher order CGPTs. The remainder 
of this section is devoted to such stability analysis. 



3.1 Analytical formula in the concentric setting 

To simplify the analysis, we assume that the receivers (sources) are evenly distributed along 
a circle of radius R centered at z. That is, 6 r = 2irr/N, r = 1,2, ... ,N, and R r = R. In 
this setting, we have A = CD, where C is an N x 2K matrix constructed from the block 
C rm = (cosm8 r smm6 r ) and D is 2K x 2K diagonal matrix: 

C 1K \ /h/R \ 



/Cn 

C21 



C12 
C22 



>2K 



\Cjvi C 



iV2 



CnkJ 



V(2i? 2 ) 



V 



I 2 /(KR«)J 



Here I2 is the 2x2 identity matrix. We note that C and D account for the angular 
and radial coefficients in the expansion of MSR, respectively. The matrix C satisfies the 
following important property; see Appendix A. 3. 

Proposition 3.1. Suppose that 2K < N holds. Then 



C*C = — \ik ■ 
2 



(3.3) 



Henceforth, we assume that the number of receivers is large enough so that 2K < N. 
In this setting, the least squares problem (3.2) admits an analytical expression as follows. 

Lemma 3.2. In the above concentric setting with sufficiently many receivers, i.e., 2K < N , 
the least squares estimation (3.2) is given by 



M 



est 



(A^d-i^vcd- 1 . 



(3.4) 



Proof. Firstly, (3.3) implies that A has full rank, so ker(L) = {0}. Moreover, 



Hence, 



which yields (3.4). □ 



M 1 



est 



(A* A)- 1 = -D- 2 . 



(— ) 2 D- 2 DC*VCDD 
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3.2 Measurement noise and stability analysis 

We develop in the rest of this section a stability analysis for the least squares reconstruction 
of CGPT from the MSR matrix, in the setting of concentric receivers (sources). 

Counting some additive measurement noise, we modify the expression of MSR to 

V = CDMDC* + E + cr noisc W. (3.5) 

Here, E is the truncation error due to the finite order K in expansion (2.11), W is an 
N x N real valued random matrix with independent and identically Gaussian entries with 
mean zero and unit variance, and <T no ise is a small positive number modeling the standard 
deviation of the noise. 

Recall that the unknown M consists of CGPTs of order up to K of the relative domain 
SB = D — z, where 5 denote the typical length scale of the domain D. The receivers and 
sources are located along a circle of radius R centered at z. Let e = 5/R be the ratio 
between the two scales, and it is assumed to be smaller than one. Due to the scaling 
property of CGPT (see (4.3)), the entries of the CGPT block M mn (8B) is 5 m+n M mn {B). 
Consequently, the size of V itself is of order e 2 , which is the order of the first term in the 
expansion (2.19). The truncation error E is of order e K+2 ; see Appendix A.l. 

According to the above analysis, we assume that the size of the noise satisfies 

Ne K+2 < a noisc < e 2 . (3.6) 

This is the regime where the measurement noise is much smaller than the signal but much 
larger than the truncation error. The presence of N in (3.6) will be clear later; see remark 
3.4. We define the signal-to-noise ratio (SNR) to be 



SNR 



c 2 



We will investigate the error made by the least squares estimation of the CGPT matrix, in 
particular the manner of its growth with respect to the order of the CGPTs. Given a SNR 
and a tolerance number tq, we can define the resolving order mo to be 



m = min <1 <m< K : J ^ ™ m F < t > . (3.7) 

y 1 1 rvJ-mm 1 1 p 

We are interested in the growth of tuq with respect to SNR. 

We have used the notation M mn , m,n = 1,...,K, to denote the building block of 
the CGPT matrix M in (2.20). In the following, we also use the notation (M)jfc, j,k = 
1, . . . , 2K, to denote the real valued entries of the CGPT matrix. 

Theorem 3.3. Assume that the condition of Lemma 3.2 holds; assume also that the additive 
noise is in the regime (3.6), Then for j,k so that (MW is non-zero, the relative error in 
its reconstructed CGPT satisfies 



/ E|(M^) jfc -( M ) jfc | 2 (Tnoiso [7 y 2l _ rfe/2l 

l(M), fc | 2 - N 



~ f 




~k~ 


2 




2 



(3. 
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Here, the symbol \l~\ is the smallest natural number larger than or equal to I. For vanishing 
(M)jfc, the error y / ¥,\(M. cst )jf s — (M)^ 2 can be bounded by the right-hand side above with 
e replaced by R~ l . In particular, the resolving order mo satisfies 

(mQE 1 "™ 10 ) 2 ~ ToSNR, (3.9) 

where tq is the tolerance number. 

Proof. From the analytical formula of the least squares reconstruction (3.4) and the 
expression of V (3.5), we see that for each fixed j, k = 1, . . . , 2K, 

(M cst - M) jk = ^(D^C'WCD- 1 ),, + ^(D-'C'ECD^ 1 )^. 

Let us denote these two terms by Ijki and Ijk2 respectively. For the first term, define W 
to be (y2/7V C)'W(y / 2/iVC), which is an N x N random matrix. Due to the orthogonality 
(3.3), W remains to have mean zero Gaussian entries with unit variance. Because D is 
diagonal, we have for each j,k = 2K, 



E(X, fcl ) 2 = ^|i-(D^)^E|W jfc | 2 (D fcfc )^ 2 = 267r4 j|° iSe ^ 2(rj/2l + rfc/21 



~3~ 


2 


~k~ 


2 




2 



Note that \j/2}\k/2] is the order of CGPT element (M) jk ; see (2.20). It is known that 
(M) jk (5B) = <yri/2l4W21(M) ifc (.B). When this ter m is non-zero, it is of order 5^/21+^/21. 
This fact and the above control of Ijki show that y / E^fc7F7T(^)j^F satisfies the estimate 
in (3.8). 

For the second term, since E is symmetric, it has the decomposition E = P*£P, where P 
is an NxN orthonormal matrix, and E is an N x N diagonal matrix consisting of eigenvalues 
of E. Then ( v / 27iVC)*E( x /2/iVC) can be written as Q*£Q where Q = t/PPC is an 
N x IK matrix satisfying Q*Q = I2K- Then the calculation for Ijki shows that 



2 6 7T 4 ' N 



Cm? = ^R^ww 



~ f 


2 


~k~ 




2 




2 





Since E is of order e K+2 as shown in (A.l), the sum is of order Ne K+2 . Therefore, we have 



E\x jk2 \* < c^+ 2 -ri/ 2 i-r^i r§l r§i - 

Since we assumed that (3.6) holds, this error is dominated by the one due to the noise. 
Hence, (3.8) is proved. 

For diagonal blocks M mm , their Frobenius norms do not vanish and (3.7) is well defined. 
In particular, (3.8) applied to the case j, k = 2m — 1, 2m, shows that the relative error made 
in the block M. mm is of order o- no \ sc m 2 e~ 2m . Using the definition of SNR, we verify (3.9). 
□ 

Remark 3.4. If E has only several (of order one) non-zero eigenvalues, then the preceding 
calculation shows that (Ijk2) 2 < Ce 2 ^ K+2 ^ and condition (3.6) can be replaced with e K+2 

^noise ^ £ • 
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(4.1) 



4 Complex CGPTs under Rigid Motions and Scaling 

As we will see later, a complex combination of CGPTs is most convenient when we con- 
sider the transforms of CGPTs under dilatation and rigid motions, i.e., shift and rotation. 
Therefore, for a double index ran, with m,n = 1, 2, . . ., we introduce the following complex 
combination of CGPTs: 

N« (A, D) = (M% n - i\CJ + i(M£ n + ilCJ, 
(A, £>) = (M- n + M-J + i(M- n - M-J. 

Then, from (2.9), we observe that 

N£(A,D)= / P n (y)(A/-/Cl ) )- 1 [(^VP m )](y) ( i S (?/), 

JdD 

N^(A,D)= / P n (y)(A/-/Cl ) )- 1 [^,V7^)](y)d S (y), 
•/an 

where P n and P m are defined by (2.12). In order to simplify the notation, we drop A in the 
following and write simply Nmn(D), N^(D). 

We consider the translation, the rotation and the dilatation of the domain D by intro- 
ducing the following notation: 

• Shift: T Z D = {x + z, for x G D}, for z G R 2 ; 

• Rotation: R e D = {e ie x, for x G £>}, for (9 G [0,2vr); 

• Scaling: sD = {sx, for x G D}, for s > 0. 

Proposition 4.1. For aZZ integers m, n, and geometric parameters 6, s, and z, the following 
holds: 

Ngl(ReD) = e^ m+n ^l(D), N$ n (R e D) = e^^N^D), (4.2) 
NW( S D) = S ™+"NW(D), N|£( Sj D) = s m+n NW(D), (4.3) 



KM(T„D) = E E C z ml NgHD)C* nk , N% n (T z D) = J2H CZ m.Mk(D)C z nk , (4.4) 



(4.5) 



1=1 fc=i i=i fc=i 

where C z is a lower triangle matrix with the m, n-th entry given by 

rn 
n 

and C z denotes its conjugate. Here, we identify z = (z±, z^) with z = z\ + iz2- 

An ingredient that we will need in the proof is the following chain rule between the 
gradient of a function and its push forward under transformation. In fact, for any diffeo- 
morphism T from R 2 to R 2 and any scalar- valued differentiable map / on R 2 , we have 



d{foT)l(h) = [df\ T{x) odTl){h), (4.6) 



10 



for any tangent vector li£l 2 , with dT being the differential of T. 

Proof of Proposition J^.l. We will follow proofs of similar relations that can be found in 
[4]. Let us first show (4.2) for the rotated domain Dg := RgD. For a function <p(y),y 6 dD, 
we define a function ip e (yg), yg := Rgy £ dDg by 

V 6 {ye) = <p° R-e(ye) = <p(y)- 

It is proved in [4] that XI — K,* D is invariant under the rotation map, that is, 

(XI - Kb a )[<p ](y g ) = (XI - KhM(v). (4.7) 

We also check that P m (R e y) = e im9 P m (y). 

We will focus on the relation for Nmn, the other one can be proved in the same way. 
By definition, we have 

N« (D) = [ P n (y)^ D , m (y)ds(y), 



dD (U 



N mn(A0 = / Pn(ye)<PD g ,m(ye)ds(y e ), 

J 8Dg 

where 

<PD, m {y) = (AI-ZCBrV, VP m )](y), 

<PD e ,m(Ve) = {M-K De r l [(v,VP m )](y e ). 

Note that the last function differs from ip s D . By the change of variables ye = Rgy in the 
first expression of (4.8), we obtain 

N« (D) = / P n (R-oy 9 ) ¥>D,m(R-6ye)ds(yg) 

J dDg 



e -ine 



p n(ye)^D,m(ve)ds(ye 

dDg 



From (4.7), we have 

(XI - JC* D M, m ](ye) = (XI - JC D )[p D , m ]{y) 
= (v y ,VP m (y)). 

Moreover, P m (y) = e~ tm6 P m (ye) so that, by applying the chain rule (4.6) with / = P m , 
T = Rg, x = y and h = u y , we can conclude that 

(v y ,VP m (y)} = e- ime (Rgv y ,VP m (yg)) 
= e- ime (v yg ,VP m (yg)). 

Therefore, (p 9 D m = e- im9 ip De , m , and we conclude that N%H(D g ) = e< m+n ^N^l(D). 

The second identity in (4.2) results from the same computation as above (the minus 
sign comes form the conjugate in the definition of N^ 2 )), and the two equations in (4.3) are 
proved in the same way, replacing the transformed function (p e by 

tp s (sy) = tp(y). 
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Thus, only (4.4) remains. Since the difference between these two comes from the conju- 
gation, we will focus only on the first identity in (4.4). The strategy will be once again the 
following: for a function <p(y),y G dD, we define a function (p z (y z ),y z = y + z € dD z , with 
D z := T Z D, by 

<P z (yz) = ( P° T- Z {y z ) = ip(y), 
which also verifies an invariance relation similar to (4.7) 

(XI - K* Dz W]{y z ) = (XI - Kh)[<p](v). (4.9) 
Moreover, for every integer q € N one has the following 

P q (y z ) = (y + z) q = E ( ff W~ r - ( 4 - 10 ) 

r=0 

Equations (4.8) become 

N« (D) = [ P n (y)<p D>m (y)ds(y), 

J 3D 

(D z ) = [ P n (yz)VD z ,m(y z )ds(y z ), 
JdD z 

where 

VdAv) = (M-IChr 1 ^, VP m )](y), 

v>D„ m (vz) = (xi-ic* Dz )- x [{^P m )]{yz). 

Thus, combining (4.9) and (4.10) leads us to 

(XI - JC* D J[(p Dzim ](y z ) = (uy z ,VP m (y z )) 



1=1 ^ ' 

Y,n)z m - l (Xi-lCh)[^i}(y) 

i=i ^ ' 

J^(^y m -\XI-Kl z W D ^(y z ), 



so that we have 



¥D z ,m(y) = E [l) 2 " 1 l( fD,M- 
1=1 ^ ' 

Hence, returning to the definition of Nmi(f 2 ) with the substitution y z o y, we obtain 

nW (d z ) = jr (j)* m - l f dD Pn(yz)^ D ,i{yz)ds(y z ), 



1=1 

m n 



1=1 k=l v / v / 



which is the desired result. Note that the index k begins with k = 1 because f dD (p z D l = 0. 
This completes the proof. □ 
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4.1 Some properties of complex CGPTs 

We define the complex CGPT matrices by := (N^) m> „ and N^ 2 ) := (N^) m , n . We 
set w = se l9 and introduce the diagonal matrix G w with the m-th diagonal entry given by 
s m e ime . Proposition 4.1 implies immediately that 

N^(T z sR e D) = C z G w N {1 \D)G w (C z y, (4.11) 

N {2 \T z sR e D) = C z G w N {2) (D)G w (C z ) t , (4.12) 

where C z is defined by (4.5). Relations (4.11) and (4.12) still hold for the truncated CGPTs 
of finite order, due to the triangular shape of the matrix C z . Using the symmetry of the 
CGPTs ([7, Theorem 4.11]) and the positivity of the GPTs as proved in [7], we easily 
establish the following result. 

Proposition 4.2. The complex CGPT matrix is symmetric: (N^ 1 ))' = N^ 1 ), and N^ 2 ) 
is Hermitian: (N^) H = N^ 2 ). Consequently, the diagonal elements of N^ 2 ) are strictly 
positive if X > and strictly negative if X < 0. 

Furthermore, the CGPTs of rotation invariant shapes have special structures: 

Proposition 4.3. Suppose that D is invariant under rotation of angle 2ir/p for some integer 
p>2, i.e., R 2lT /pD = D, then 

N^(D) = 0, if p does not divide (m + n), (4-13) 

N$ n (D) = 0, ifp does not divide (m - n). (4-14) 

Proof. Suppose that p does not divide (m + n), and define r := 2ir(n + m)/p mod 2tt. 
Then by the rotation symmetry of D and the symmetry property of the CGPTs, we have 

NW (D) = N« (R 2w/p D) = e« m+n ^N£ n (D) = e *N« (JD). 
Since r < 2ir and r / 0, we conclude that Nmi(D) = 0. The proof of (4.14) is similar. □ 

5 Shape Identification by the CGPTs 

We call a dictionary V a collection of standard shapes, which are centered at the origin 
and with characteristic sizes of order 1. Given the CGPTs of an unknown shape D, and 
assuming that D is obtained from a certain element B € T> by applying some unknown 
rotation 9, scaling s and translation z, i.e., D = T z sRgB, our objective is to recognize 
B from T>. For doing so, one may proceed by first reconstructing the shape D using its 
CGPTs through some optimization procedures as proposed in [11], and then match the 
reconstructed shape with T>. However, such a method may be time-consuming and the 
recognition efficiency depends on the shape reconstruction algorithm. 

We propose in subsections 5.1 and 5.2 two shape identification algorithms using the 
CGPTs. The first one matches the CGPTs of data with that of the dictionary element by 
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estimating the transform parameters, while the second one is based on a transform invariant 
shape descriptor obtained from the CGPTs. The second approach is computationally more 
efficient. Both of them operate directly in the data domain which consists of CGPTs and 
avoid the need for reconstructing the shape D. The heart of our approach is some basic 
algebraic equations between the CGPTs of D and B that can be deduced easily from (4.11) 
and (4.12). Particularly, the first four equations read: 

N { $(D)=w 2 N$(B), (5.1) 
N$(Z>) = 2N[ 1 »(D)z + w 3 N®(B), (5.2) 
N®(D)=*N®(B), (5.3) 
IN© (£>) = 2N^>(D)z + s 2 wN ( $(B), (5.4) 

where w = se lS . 

5.1 CGPTs matching 

5.1.1 Determination of transform parameters 

Suppose that the complex CGPT matrices NW(5),N( 2 I(B) of the true shape B are given. 
Then, from (5.3), we obtain that 



Ng>(D)/Ng>(JS). (5.5) 
Case 1: Rotational symmetric shape. If the shape B has rotational symmetry, i.e., 

(2) 

R 27T / p B = B for some p > 2, then from Proposition 4.3 we have N^ 2 (B) = and the 
translation parameter z is uniquely determined from (5.4) by 

NfJ(D) , . 

(5.6) 



2N< 2 1 ) (£>)' 



On the contrary, the rotation parameter 9 (or e %e ) can only be determined up to a multiple 
of 2tt/p, from CGPTs of order [p/2] at least. Although explicit expressions of e ipe can be 
deduced from (5.1) - (5.4) (or higher order equations if necessary), we propose to recover 
e ip9 Dy so iving the least squares problem: 



mm 



(jN^(T z sR e B) - NW(D)||^ + ^ 2 \T z sR e B) - N( 2 )(£>)|Q . (5.7) 

Here, s and z are given by (5.5) and (5.6) respectively, and 'N^(D) and ~N( 2 \D) are the 
truncated complex CGPTs matrices of dimension [p/2] x [p/2]. 

Case 2: Non rotational symmetric shape. Consider a non rotational symmetric shape 
B which satisfies the assumption: 

N^fij^O and det( N |i! (B) N || (jB) j + 0. (5.8) 
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From (5.2) and (5.4), it follows that we can uniquely determine the translation z and the 
rotation parameter w = e td from CGPTs of orders one and two by solving the following 
linear system: 

Ng(D)/Ng>p) = 2z + wN^(B)/N^(B), 

N%>{D)/N<8(D) = 2z + W N^{B)/N^(B). (5.9) 



5.1.2 Debiasing by least squares solutions 

In practice (for both the rotational symmetric and non rotational symmetric cases), the 
value of the parameters z, s and provided by the analytical formulas and numerical proce- 
dures above may be inexact, due to the noise in the data and the ill-conditioned character 
of the linear system (5.9). Let z*,s*,9* be the true transform parameters, which can be 
considered as perturbations around the estimations z, s, obtained above: 



z + 5 Z , s* = s5 s , and 6* = 6 + 5$ 



(5.10) 



for 5 z ,5g small and 5 S close to 1. To find these perturbations, we solve a nonlinear least 
squares problem: 



mm 

*i «/ a' 



N {1 \T z is'R d >B) - N (1) (D) + N {2) (T z , s' R g/ B) - N (2) (£>) 



(5.11) 



with (z,s,9) as an initial guess. Here, the order of the CGPTs in (5.11) is taken to be 
2 in the non rotational case and max(2, [p/2]) in the rotational symmetric case. Thanks 
to the relations (4.11) and (4.12), one can calculate explicitly the derivatives of the objec- 
tive function, therefore can solve (5.11) by means of standard gradient-based optimization 
methods. 



5.1.3 First algorithm for shape identification 

For each dictionary element, we determine the transform parameters as above, then measure 
the similarity of the complex CGPT matrices using the Frobenius norm, and choose the 
most similar element as the identified shape. Intuitively, the true dictionary element will 
give the correct transform parameters hence the most similar CGPTs. This procedure is 
described in Algorithm 1. 

5.2 Transform invariant shape descriptors 

From (5.3) and (5.4) we deduce the following identity: 

NSW + ,5,2, 



2Nf 1 »(B) 2NS?(fl) 

which is well defined since ^ thanks to the Proposition 4.2. Identity (5.12) shows a 
very simple relationship between — and — ^y for D = T z sRgB. . 

2N 1:L (B) 2N 1:L (D) 
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Algorithm 1 Shape identification based on CGPT matching 

Input: the first A;-th order CGPTs N^^D^N^^) of an unknown shape D 
for B n £ V do 

1. Estimation of z,s,9 using the procedures described in subsections 5.1.1 and 5.1.2; 

2. D <- R_ e s- l T^ z D, and calculate N^(D) and N( 2 )(l)); 

3. EM <- N( 1 )(S n ) - N^(D), and E^ <- N^(B n ) - N^(D); 

4. e n (\\EW\\ 2 F+ ||^ 2 )||^) 1 / 2 /(||N( 1 )(i? n )||^+ \\N^(B n )\\ 2 F )^; 

5. n <— n + 1; 
end for 

Output: the true dictionary element re* <— argmin n e n . 



Let u = — m We first define the following quantities which are translation invariant: 

2Nf 1 ) (Z)) bH 

J^\D) = N«(T_ U D) = C- U NW(D)(C- U )*, (5.13) 
j( 2 )(D) = N (2) (T_ U D) = CF«N (2) (D)(Cr")*, (5.14) 

with the matrix C~ u being the same as in Proposition 4.1. From J^(D) = {Jmrn{D)) m ^ n 

and j( 2 \D) = (Jmm(D)) mi „, we define, for any indices m, n, the scaling invariant quanti- 
ties: 

S ^{D) = - - J - iD ? /2 , SH(D) = - (5.15) 



Finally, we introduce the CGPT-based shape descriptors = (Imn)m n and 1^ = 
(1 {2) ) ■ 

l£l(D) = \SW(D)\, l£l(D) = \S£1(D)\, (5.16) 

where | • | denotes the modulus of a complex number. Constructed in this way, andX( 2 ) 
are clearly invariant under translation, rotation, and scaling. 

It is worth emphasizing the symmetry property, Zmn = Inm,Imn = Znm, and the fact 
that Xmli = 1 for any m. 



5.2.1 Second algorithm for shape identification 

Thanks to the transform invariance of the new shape descriptors, there is no need now for 
calculating the transform parameters, and the similarity between a dictionary element and 
the unknown shape can be directly measured from 1^ and X^ 2 \ As in Algorithm 1, we 
use the Frobenius norm as the distance between two shape descriptors and compare with 
all the elements of the dictionary. We propose a simplified method for shape identification, 
as described in Algorithm 2. 
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Algorithm 2 Shape identification based on transform invariant descriptors 

Input: the first k-th. order shape descriptors X^\D) ,X^ 2 \D) of an unknown shape D 
for B n e V do 

1. e n ,^ (IIXW^-ZW^II^+II^ 2 )^)-^ 2 )^)^) 172 ; 

2. n <— n + 1; 
end for 

Output: the true dictionary element n* argmin n e n . 



6 Numerical Experiments 

In this section we present a variety of numerical results on the theoretical framework dis- 
cussed in this paper in the context of target identification from noisy MSR measurements. 
Given a shape Dq of characteristic size 5, the procedure of our numerical experiment can 
be summarized as follows: 

1. Data simulation. N sources (and also receivers) are equally distributed on a circle of 
radius R, which is centered at an arbitrary point zq £ -Do an d includes Dq, see Figure 
1. The MSR matrix is obtained by evaluating numerically its integral expression (2.7) 
then adding a white noise of variance cr 2 oise . For simplicity, here we suppose that the 
reference point zq £ Dq can be estimated by means of algorithms such as MUSIC 
(standing for Multiple Signal Classification) [2, 7]. 

2. Reconstruction of the CGPTs of D = Dq — zq using formula (3.4) or the least squares 
algorithm (3.2). 

3. For a given dictionary T>, apply Algorithm 1 (or Algorithm 2) using the CGPTs of D 
and identify the true shape from T>. 

We emphasize that the reconstructed CGPTs of shape D depend on the reference point zq. 
We fix the conductivity parameter k = 4/3 throughout this section. 



6.1 Reconstruction of CGPTs 

The theoretical analysis presented in section 3 suggests the following two step method for 
the reconstruction of CGPTs. First we apply (3.4) (or equivalently solve the least squares 
problem (3.2)) by fixing the truncation order K as in (3.6): 

^<minf l0g((Tn0ise/JV) -2,iV/2V (6.1) 

V log £ J 

Here, cr no i S c is the standard deviation of the measurement noise and e = 5/R with 5 being 
the characteristic size of the target and R the distance between the target center and the 
circular array of transmitters/receivers. Then, we keep only the first rriQ orders in the 
reconstructed CGPTs, with rriQ being the resolving order deduced from estimation (3.9): 

log Cr noisc - log T 

m o = ^ , (6.2) 

21oge 
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Figure 1: An example of the configuration for MSR data simulation. The unknown shape 
is an ellipse whose long and short axes are 2 and 1, respectively. N = 51 sources/receivers 
(marked by "x") are equally placed on a circle of radius R = 2 centered at zq = [0,0] 
(marked by "*")■ 

and ro < 1 is the tolerance number introduced in (3.7). In all our numerical experiments 
we set the noise level <r no ise 

to: 

"noise — (Vmax - VrmnKo) (6.3) 

with a positive constant <jq and V max and V m j n being the maximal and the minimal co- 
efficient in the MSR matrix V. Using the configuration given in Figure 1 and for various 
noise level, we reconstruct the CGPTs of the ellipse up to a truncation order K which is de- 
termined as in (6.1). For each k < K, the relative error of the first k-th order reconstructed 
CGPTs is evaluated by comparing with their theoretical value ([7, Proposition 4.7]). The 
results are shown in Figure 2. In Figure 3 we plot the resolving order m,Q given by (6.2) 
and the relative error of the reconstruction within this order, for do in the range [10 -3 , 1]. 

6.2 Dictionary matching 

We are now ready to present the results of the dictionary matching algorithms discussed in 
the sections 5.1 and 5.2. Unless specified, in the following we suppose that the unknown 
shape of the target Do is an exact copy of some element from the dictionary, up to a rigid 
transform and dilatation. As examples, we consider a dictionary of flowers and a dictionary 
of Roman letters. The aim is to identify the target D$ from imaging data if it belongs to 
one of the dictionaries. 

6.2.1 Matching on a dictionary of flowers 

We start by considering a simple dictionary of rotation invariant "flowers", on which the 
shape identification algorithm can be greatly simplified. The boundary of the p-th flower 
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10 1 1 ' 

2 4 6 8 10 12 14 16 18 5 10 15 

Order of CGPT Order of CGPT 



(a) cro = 0.01, mo = 6 (b) ctq = 0.1, m = 4 




Order of CGPT Order of CGPT 

(c) cro = 0.5, mo = 3 (d) era = 1.0, mo = 2 

Figure 2: Relative error of the reconstructed CGPTs. For each noise level, we repeat the 
experiment 100 times (corresponding to 100 realizations of the noise) and the reconstruction 
is taken as their mean value. The horizontal solid line in each figure indicates the resolving 
order mo given by (6.2) with the tolerance number ro = 10 _1 . 
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10" 3 10~ 2 10"' 10° 10' 3 10" 2 10"' 10° 



sigma sigma Q 

(a) Resolving order (b) Relative error 

Figure 3: The resolving order mo, for <to £ [10~ 3 , 1],to = 10 _1 , and the relative error of the 
reconstruction within this order. As in Figure 2, we repeat the experiment 100 times and 
the reconstruction is taken as their mean value. The large variations of the relative error 
in (b) for <To > lO^ 1 indicate the instability of the reconstruction for very noisy data. 



Bp is defined as a small perturbation of the standard disk: 

dB p (0 = x(f)(l +Vcos(pO), x® = , (6.4) 

where p > 2 is the number of petals and r] > is a small constant. According to Proposition 
4.3, Nmn(Bp) is zero if p does not divide m + n. For an unknown shape D = T z sRqB p , the 

translation parameter is given by z = ^(2) ^ ■ Moreover, simple calculations show that 

ZW(D) and ~N^(B p ) have exactly the same zero patterns. 

Therefore, we can find the true number of petals by searching the first nonzero anti- 
diagonal entry in I^(D). 

We fix T] = 0.3 (the amplitude of the perturbation introduced in (6.4)) and 5/R = 
0.5. The unknown shape Dq is obtained by applying the transform parameters z = 
[16.3, —46.7], s = 7.5,9 = 2.69 on B p , and the reference point for data acquisition is 
zq = [15, —45.5]. The results for two flowers of 5 and 7 petals are shown in Figure 4, where 
we plot the mean absolute value of the anti-diagonal entries mn, for m+ n = I, I = 2, . . . , 11, 
in (D) by varying the noise level a®. One can clearly distinguish the peak which indicates 
the true number of petals for o"o up to 10~ 2 . 

Stability. Let us consider now the model (6.4) with a general C 1 function h(£) in place 
of cos(p£). It was proven in [4] that: 

N^ n (B p ) = 27r ?? ^X +n + 0( V 2 ). (6.5) 

Therefore as long as the perturbation h(£) is close to cos(p£), the significant nonzero coeffi- 
cients in I^(D) will concentrate on the same anti-diagonals. We confirm this observation 
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23456789 10 11 23456789 10 11 



(a) p = 5 (b) p = 7 

Figure 4: Mean values of the anti-diagonal entries of for the flowers of 5 and 7 petals 
at different noise levels. 



by applying the same procedure above on a flower with one damaged petal: 

f for£e [0,2vr/p), 

8B p {0 = { (6.6) 
[x(£)(l + »7cos(p£)) for f G [27r/p,27r). 

Here, f(-,t) : M i— )• M is a polynomial of order 6, constructed such that <9-B p is C 2 -smooth, 
and t € (0, 1) is the percentage of the damage; see Figure 5. In Figure 6 we plot the mean 
value of the anti-diagonal entries at different noise levels. Compared to Figure 4, we see 
that the effect of the damage in the petal dominates the measurement noise. Nonetheless, 
the peak indicating the true number of petals is still visible. 




(a) r = 0.5 (b) r = 0.8 

Figure 5: Flowers with one damaged petal. The following parameters are used in (6.6): 
p = 7 iV = 0.3, t = 0.5 for (a) and t = 0.8 for (b). 
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n+m n+m 



(a) To = 0.5 (b) r = 0.8 

Figure 6: Mean value of the anti-diagonal entries of Z^ 1 ) for the flowers of Figure 5 at 
different noise levels. The peaks indicate the number of petals. 

6.2.2 Dictionary of letters 

Next we consider here a dictionary consisting of 26 Roman capital letters without rotational 
symmetry. The shapes are defined in such a way that the holes inside the letters are filled, 
see Figure 11. We set 5/R = 0.5, s = 2.4762,6* = 6.0827,2 = [33.3505,73.8395] and the 
center of mass of the target at [33.4042, 73.8627]. 

Performance of Algorithm 1. First we test Algorithm 1 on the letter "P". For the 
noiseless case (<ro = 0), the values of e„ defined in Algorithm 1 are plotted in Figure 7 
(a) and (b). These results suggest that the high order CGPTs can better distinguish 
similar shapes such as "P" and "R", since they contain more high frequency information 
[4]. Nonetheless, the advantage of using high order CGPTs drops quickly when the data are 
contaminated by noise, and low order CGPTs provide more stable results in this situation, 
see Figure 7 (c) and (d). 

By repeating the same procedure as above, we apply Algorithm 1 on all letters at noise 
levels do = and ctq = 0.1, and show the result in Figure 8 (a) and (c). At the coordinate 
(m, n), the unknown shape is the m-th letter and the color represents the relative error 
(in logarithmic scale) of the CGPTs when compared with the ra-th standard letter of the 
dictionary. 

Stability. In real world applications we would like to have Algorithm 1 work also on 
letters which are not exact copies of the dictionary, such as handwriting letters. Figure 12 
shows the letters obtained by perturbing and smoothing the dictionary elements. With 
these letters as unknown shapes, we repeat the experiment of Figure 8 (a) and (c) by 
applying Algorithm 1 on the standard dictionary and show the results in Figure 8 (b) and 
(d). Comparing with the results of Figure 8 (a) and (c), we see that Algorithm 1 remains 
quite stable, despite of some slight degradations. 
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(a) o"o = 0, order < 2 



(c) o-o = 0.1, order < 2 





(b) (Jo = 0, order < 5 



lllllll.h 



.11.1 



(d) (To = 0.1, order < 5 



Figure 7: The identification of the letter "P" using the first 2, and 5 orders CGPTs at noise 
levels o"o = and do = 0.1. The bar represents the relative error e n between the CGPTs 
of the n-th letter and that of the data, as defined in Algorithm 1, and the shortest one in 
each figure corresponds to the identified letter. For (c) and (d), the experiment has been 
repeated for 100 times, using independent draws of white noise, and the results are the 
mean values of all experiments. 
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Dictionary letter 

(a) <7o = 0, order < 5, Standard letters 



Dictionary letter 

(b) (Jq = 0, order < 5, Perturbed letters 




Dictionary letter 

(c) (To =0.1, order = 1, Standard letters 



Dictionary letter 

(d) (Jo = 0.1, order = 1, Perturbed letters 



Figure 8: Algorithm 1 applied on the all 26 letters using the standard dictionary (Figure 11) 
at noise level ctq = (first column) and <ro = 0.1 (second column), with the color indicating 
the relative error e n in logarithmic scale. The unknown shapes in the first row are exact 
copies of the standard dictionary, and in the second row are those of Figure 12. In (a) 
all letters are correctly identified, while in (b) letters 'E' is identified as 'H'. For the noisy 
case, the experiment has been repeated 100 times, using independent draws of white noise, 
and the results in (c) and (d) are the mean values of all experiments, where only the first 
order CGPT is taken into account. 22 and 21 letters are correctly identified in (c) and (d), 
respectively. 
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Performance of Algorithm 2. In the case of noiseless data, Algorithm 2 provides correct 
results with low computational cost. Here we repeat the experiment in Figure 7 (a) and (c) 
using Algorithm 2, and plot the error e n defined in Algorithm 2 in Figure 9. Nonetheless, 
when data are noisy, Algorithm 1 performs significantly better than Algorithm 2, as shown 
by Figure 10 where we compare the two algorithms for identifying letter "P" at various 
noise levels. Thanks to the debiasing step (5.11), Algorithm 1 is much more robust with 
respect to noise than Algorithm 2, in which there is no debiasing and the invariance of the 
shape descriptors and 1^ may be severely affected by noise (see Figure 10). 
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Dictionary letter 

(a) o"o = 0, order < 5, Standard letters 



Dictionary letter 

(b) <7o = 0, order < 5, Perturbed letters 



Figure 9: Algorithm 2 applied on the all 26 letters using the standard dictionary (Figure 11) 
at noise level <7o = 0. The unknown shapes in (a) are exact copies of the standard dictionary, 
while in (b) are those of Figure 12. The color indicates the error e n in logarithmic scale. 
All letters are correctly identified in both (a) and (b). 



7 Conclusion 

In this paper, we have designed two fast algorithms which identify a target using a dictionary 
of precomputed GPTs data. The target GPTs are computed from multistatic measurements 
by solving a linear system. The first algorithm matches the computed GPTs to precomputed 
ones (the dictionary elements) by finding rotation, scaling, and translation parameters and 
therefore, identifies the true target shape. The second algorithm is based on new invari- 
ants for the CGPTs. We have provided new shape descriptors which are invariant under 
translation, rotation, and scaling. The stability (in the presence of additive noise in multi- 
static measurements) and the resolution issues for both algorithms have been numerically 
investigated. The second algorithm is computationally much cheaper than the first one. 
However, it is more sensitive to measurement noise in the imaging data. To the best of our 
knowledge, our procedure is the first approach for real-time target identification in imaging 
using dictionary matching. It shows that GPT-based representations are an appropriate 
and natural tool for imaging. Our approach can be extended to electromagnetic and elastic 
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Figure 10: Comparison of Algorithm 2 and Algorithm 1 on identification of the standard 
letter "P". At each noise level, the experiment has been repeated 1000 times, using inde- 
pendent draws of white noise. For each algorithm, the curve represents the percentage of 
experiments where the letter "P" is correctly identified. 

imaging as well [12, 5]. We also to plan to use it for target tracking from imaging data. 



A Appendix: Several Technical Estimates 
A.l The truncation error in the MSR expansion 

Recall the expansion of the element in the MSR matrix (2.11). We prove the following 
estimate of the truncation error. 

Proposition A.l. Let E rs be as in (2.11). Set e = 5/R, the ratio between the typical 
length scale of the inclusion D and the distance of the receivers (sources) from the inclusion. 
Assume also that e is much smaller than one. Then 



\E r A<e K+ \ 



(A.l) 



Proof. From the Taylor expansion of multivariate functions ([24], Chapter 1), we verify 
that the truncation error E rs can be written as 



/ e K (y,x r ,z)(XI - K* D ) 1 

JdD 



dT(- - x a ) 
dv 



(y)ds(y) 



[ r x (y;x r ,z)(A/-/Cl ) )- 1 

JdD 



d_ 

dv 



e K (",z,x s ) 



(y)ds(y). 
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Here, Txiy] x r , z) and e^(y; x r , z) (and similarly e^(y; z, x s )) are given by 

(-l)l' 



T K (y;x r ,z)=Y J E ^r^^r(x r - z)(y - zf, 

fc=l |a|=fc 



|a|=i<-+l ^° 

Due to the invariance relation (4.7), the operator (XI — /C^) _1 , as an operator from the 
space L 2 (dD) to itself, is bounded uniformly with respect to the scaling of D. Consequently, 
the first term in E rs is bounded by 

dT(- — x s ) i dT(- — x s ) 
C\\eK{-,x r , z)\\ L oo( dD )\\ — \\L 2 (dD)\9D\2 < C \\eK\\L°°(dD)\\ \\L°°{dD)\dD\. 

Assume that z £ D; the distance between D and the receivers (sources) is of order R. From 
the above expression of ex, the explicit form of d a T in (2.18), and the fact that \y — z\ < C5 
for y £ D, we have 



Mj/;s r ,z)| <C ( £ l.\\d a T r (x r 

\\a\=K+l a ' 



\c(d) i \v-A K+l ^ C \r 



Similarly, we have 1 1 c?^ r ( - — £,s)||z,°o(,9D) < CR~ l . The measure \dD\ in dimension two is of 
order 5. Substituting these estimates into the bound for the first term in E rs , we see that 
it is bounded by Ce K+2 . 

The second term can be bounded from above by 

C || Fk I \\L°°(dD)\dD\. 

We have ||r^(-; x r , z) \\l°°(8D) < Ce, which is the order of the leading term. Further, from 
the explicit form of ex, we verify that 

\\ deK{ '^ Xs) h^8 D) < c (||r(- - x s )\\ CK+2(B) 5 K ^ + ||rc- - x s )\\ cK+1(B) 5 K ) < Cj?L. 

As a result, the above upper bound for the second term in E rs is of order £ K+2 as well. 
This proves (A.l). □ 

Proposition A. 2. The solution u s (x) defined by the transmission problem (2.2) satisfies 
the symmetry property 

u s (x r ) = u r (x s ). (A.2) 

Proof. Let Sl £ s be the the ball of radius e centered at x s , and the ball of radius e centered 
at x r . Let Q e be the domain Bn\(Q £ U Q% U D) where Br is a sufficiently large ball with 
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radius R. Then we have 



u s (x)Au r (x) — u r (x)Au s (x) )dx 



u s (x)-^-(x) - u r (x)-^-(x)Jds(x) - 



u s (x)-^-{x) - u r (x)-^-(x) )ds{x) 



am 



u s (x)-^ L (x) - u r (x)-^-(x) )ds(x) 



u s (x)— L (x) -u r (x)—^(x) )ds(x) + 
on + on 



1 3D 

Jg + + Jd + Jr- 



dB R 



u s (x) (x) - U r (x) ^ (x) 
on + on 



ds(x) 



For J/), thanks to the jump conditions in (2.2), we have that 
J D = K 



8D 



c)^^-(x) — u s (x)^^(x) jds(x) = k J ( 



The other two terms Jf and Jjf can be treated similarly; hence we focus on the first 
item. We've shown that u s (x) = T(x — x s ) + 5d[(/> s ]. In a neighborhood of £l e s , we have 



\ur\\L°° + ||Vu r || L oo + ||5d[0 s ]||l=° + ||V5d[0 s ]||l°° < C. 



Consequently, 



< C 



u r {x) 



mi 



an 

( 7^0*0 - ^-(x-x s ) 
\ on on 



dB e {x s ) 

ds(x) < 



(1 + \ \oge\)ds(x) < Ce|loge|. 



u r (x) ^^ (x)ds(x] 
dn<i 9n 



< Ce. 



These estimates imply that 



lim Jf = lim 

£->0 £-»0 



dT 



dB e (x s ) 



u r (x s + y) — {y)ds{y) = lim 



dn 



1 f 2n 

— / eu r (x s + e9)d9 = u r (x s ). 
Jo 



The same analysis applied to shows that lim £ _>.o Jr = —u s (x r )- 

To control Jr, we recall the fact that <Sz> [<?!>] decays as and VSp [<f>] decays as \x\~ 2 
for <p G L 2 (dD) satisfying f QD cj)ds = 0; these estimates imply that the logarithmic part of 
u s dominates. Therefore, 



lim Jr = lim 

R-^roo R— >oo 



1 I I ^ 1 *^ *^ T ^ i I I ^ 5 S / 7 / \ 

log |x — x s \— -= log \x — x r \— ds(x). 



The integrand above can be written as 



■ X X g \ ( t' - 3j X<yj .. I I 

lo S TZ — I Z~l2~ + lo S l x ~~ x r\ 



(v x ,x x r ) (u x ,x x s ) 



We verify that the first term is of order o(-jjj); its contribution to the limiting integral is 
hence negligible. The second term in the integrand can be further written as 



log \x — x r \ 



+ 



iy x ^x x T {x x ' sj) 
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From 



1 



1 



rp I 2 I rp 1 2 I O / nr* rp \ 



x — X 



2 



rp rp I 2 I np rp 

iXj iXj ip \ tXj tXj g 



s 



we verify that the second term in the integrand is of order (log R/R 2 ); hence its contri- 
bution to the limiting integral is also zero. To summarize, we have limft_>. 00 Jr = 0. 

From the above analysis, we take the limit e — > 0, R —¥ oo on the equality = Jf + + 
Jd + Jr and conclude that (A. 2) holds. □ 

A. 2 Proof of formula (2.18) 

Formula (2.18) is well-known. We include a proof for reader's sake. 

In order to prove (2.18), we need to find the derivative of the function log \x\. To this 
end, we consider the Taylor expansion of the logarithmic function around the point x. 
The most convenient method for this expansion is to view the space variables as complex 
numbers. For a small perturbation z of the point x (x, z G C), we calculate 



To expand the first item on the right-hand side of the above equality, we write it as log(l— — ), 
and since |^| < 1 we obtain the expansion 



Taking the conjugate, we obtain the expansion for log(x — z) — log x. Consequently, we have 



In the last equality, we understood the variable z as real variable and used the representation 
(2.13). Compare the last term of the above formula with the (real- variable) multivariate 
expansion of log \ x — z\ — log we observe that 



log \x — z\ — log \x 



— ([log(x — z) — log a;] + [log(x — z) — logx]) . 







For each double index a, we get (2.18) 
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A. 3 Proof of formula (3.3) 

The proof is a straightforward computation. The elements of the matrix C*C correspond to 
inner products of columns of the matrix C, that is, the inner products of vectors formed by 
evaluating sin and cos functions at (fci#i, . . . , kiOjy) and at (&2#l> ■ ■ ■ > ^On), where k±, k 2 = 
1, 2, . . . , K, ki + k 2 < 2K < N, and 6j = 2-Kj/N, j = 1, 2, . . . , N. When two cos vectors are 
chosen, the inner product becomes 

N 1 N 

E 1 ^-^ / , 27r(fc 1 +fc 2 )j . 2ir(k 1 +k 2 )j j . 27r(fc 1 -fc 2 )j j . 2vr(fc 1 -fc 2 )j 

cos k\6j cos fe 2 ^j = ^2^V e N + e N + e N +e ^ 

Since k\ + fc 2 is an integer less than N, the first two sums always vanish because 



\ e n = = 

, 2^(fc 1+ fc 2 ) 



1 - e l ^^r 



When k\ = k 2 , the last two sums contribute and the overall result is N/2. When k\ ^ k 2 , 
the inner products under estimation is zero according to the above observation. 

The case of inner product with sin and sin or cos and cos vectors can be similarly 
analyzed, and it can be easily seen that (3.3) holds. 
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